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Advanced  Hybrid  Modeling  of  Hall  Thruster  Plumes 


Tyler  D.  Huismann  and  lain  D.  Boyd' 

University  of  Michigan,  Ann  Arbor,  MI  48109,  USA 

Abstract:  A  hybrid  particle-fluid  method  is  applied  to  model  the  plume  from  a  6  kW 
Hall  thruster  operated  in  the  Large  Vacuum  Test  Facility  at  the  University  of  Michigan. 
The  approach  utilizes  the  direct  simulation  Monte  Carlo  method  and  the  Particle-in-Cell 
method  to  simulate  the  collision  and  plasma  dynamics  of  xenon  neutrals  and  ions.  The 
electrons  are  modeled  as  a  fluid  using  conservation  equations.  A  second  code  is  employed  to 
model  discharge  chamber  behavior  to  provide  improved  input  conditions  at  the  thruster  exit 
for  the  plume  simulation.  Simulation  accuracy  is  assessed  using  experimental  data 
previously  recorded. 

Nomenclature 

A  =  area,  m2 

Ci  =  axial  velocity  component  of  VDF,  m  s"1 

C,  =  ionization  rate  coefficient,  m3  s"1 

D  =  mean  thruster  diameter 

E  =  electric  field,  V  m'1 

e  =  electron  charge,  1.6  x  10"19  C 

g  =  magnitude  of  relative  velocity,  m  s’1 

Id  =  discharge  current,  A 

Isp  =  specific  impulse,  s 

j  =  current  density,  A  m"2 

k  =  Boltzmann  constant,  1.38  x  10"23  J  K'1 

m  =  mass,  kg 

N  =  number  of  particles 

n  =  number  density,  m"3 

p  =  pressure,  Pa 

T  =  temperature,  eV 

u  =  axial  velocity,  m  s"1 

v  =  velocity,  m  s"1 

Vd  =  discharge  voltage,  V 

et  =  ionization  energy,  eV 

k  =  thermal  conductivity,  W  (K  m)'1 

v  =  collision  frequency,  s"1 

co  =  viscosity  temperature  exponent 

<[>  =  plasma  potential,  V 

xp  =  electron  velocity  stream  function,  m'2  s"1 

Gel  =  collision  cross  section,  m2 

a-  =  electrical  conductivity,  A  (V  m)"1 

Subscripts 

a  =  anode 
d  =  discharge 
e  =  electron 
i  =  ion 
TE  =  thruster  exit 
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I.  Introduction 

HALL  thrusters  are  an  efficient  propulsion  option  for  spacecraft,  with  high  specific  impulses  making  them 
particularly  well  suited  for  low  thrust  missions.  A  primary  concern  regarding  the  use  of  Hall  thrusters  is  the  effect 
of  their  plumes  on  spacecraft  integration.  Possible  spacecraft  contamination  and  communications  interference 
emphasize  the  importance  of  accurately  analyzing  Hall  thruster  plumes.  Two  aspects  of  plume  analysis  include 
numerical  simulation  and  ground-based  experiments  conducted  in  vacuum  chambers.  One  such  facility  is  the  Large 
Vacuum  Test  Facility  (LVTF)  at  the  University  of  Michigan’s  Plasmadynamics  and  Electric  Propulsion  Laboratory 
(PEPL;  see  Figure  1).  It  is  the  focus  of  the  present  work  to  assess  plume  simulation  accuracy  through  comparisons 
of  simulation  results  with  experimental  measurements  conducted  previously  at  PEPL. 

Hall  thruster  plume  modeling  has  been  reviewed  by  Boyd1  where  it  was  determined  that  hybrid  methods  are  the 
most  successful.  In  general,  the  plume  of  a  Hall  thruster  consists  of  neutrals,  energetic  ions,  and  electrons.  Plume 
behavior  is  complicated  by  the  multiple  types  of  collisions  involved,  such  as  collisions  due  to  thermal  velocity  and 
collisions  between  neutrals  and  ions  with  charge  exchange,  as  well  as  different  physical  phenomena,  such  as  self- 
consistent  electric  fields.  Furthermore,  in  ground-based  facilities  the  presence  of  a  background  gas  must  be 
accounted  for.  Computational  methods  are  well  suited  for  plume  analysis  since  different  physical  models  can  be 
interchanged  to  allow  for  varying  degrees  of  fidelity. 

In  this  paper,  a  discharge  chamber  code,  HPHall,  is  applied  to  model  a  6  kW  Hall  thruster  in  order  to  accurately 
determine  thruster  exit  conditions.  The  HPHall  predictions  are  compared  to  semi-empirical  methods  that  were 
previously  utilized  to  estimate  thruster  exit  conditions.2  A  hybrid  simulation  method  is  then  applied  in  order  to 
investigate  the  thruster  plume  in  the  LVTF.  A  direct  simulation  Monte  Carlo  (DSMC)  method3  is  used  to  model 
collision  dynamics,  and  a  Particle-in-Cell  (PIC)  method4  is  used  to  capture  electric  field  effects.  A  detailed  fluid- 
electron  model5  is  incorporated  in  place  of  the  standard  Boltzmann  relation  for  modeling  electrons  and  is  discussed 
below.  Comparisons  between  the  thruster  exit  conditions  are  made  to  determine  effects  on  the  calculation  of  plasma 
parameters.  These  comparisons  are  also  extended  to  experimental  measurements  previously  recorded  in  the  plume, 
in  order  to  determine  the  effectiveness  of  the  various  methods  utilized. 

II.  Facility  and  Thruster 

The  LVTF  is  a  6  m  diameter,  9  m  long  cylindrical  stainless  steel  vacuum  chamber.  The  facility  is  shown  in 
Figure  1.  The  chamber  is  pumped  by  seven  CVI  TM-1200  re-entrant  cryo-pumps  which  are  surrounded  by  liquid 
nitrogen  (LN2)  cooled  shrouds.  The  chamber  has  a  total  pumping  rate  of  about  240,000  L/s  for  xenon.  This  resulted 
in  a  back-pressure  of  approximately  1.3  x  10"5  torr  for  nominal  operating  conditions.  The  current  density 
measurements  utilized  in  this  study  were  obtained  with  Faraday  cup  probes.  The  near-field  current  density  contours 
consist  of  over  64,000  individual  measurements  of  local  current  density.  For  further  details,  see  Ref.  6. 

The  thruster  tested  in  this  facility  was  a  6  kW  Hall  thruster.  The  thruster  was  operated  in  the  LVTF  over  a  range 
of  operating  modes,  with  the  nominal  condition  being  the  focus  of  this  work.  See  Table  1  for  the  specific  operating 
conditions  of  this  mode.  The  thruster  is  a  nominal  6  kW  input  power  thruster  designed  to  produce  100-600  mN  of 
thrust  at  1000-3000  s  of  Isp,  all  stated  in  Ref.  6.  The  thruster  was  equipped  with  a  LaB6  center-mounted  hollow 
cathode  that  was  operated  at  around  7%  of  the  anode  mass  flow  rate.  For  further  details  of  the  thruster  operation, 
see  the  same  reference. 


Operating  Condition 

rha 

Id 

vd 

Nominal 

20  mg/s 

20  A 

300  V 

Table  1:  6k W  Hall  thruster  operating  conditions 
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Figure  1:  LVTF  chamber 
III.  Numerical  Methods 

III.l  HPHall  Methodology 
A.  Overview 

For  the  past  decade,  the  Hall  thruster  modeling  package  HPHall  has  been  the  focus  of  much  research  and 
analysis.  HPHall  models  the  plasma  within  the  thruster  discharge  chamber  and  near-field  plume,  employing  both 
fluid  and  particle-in-cell  (PIC)  numerical  methods  on  an  axisymmetric  grid.7  The  code  has  been  found  to  be 
effective  in  creating  either  time -averaged  outputs  of  performance  data  which  could  be  utilized,  for  example,  in  an 
erosion  mode,8  or  to  examine  the  time-varying  nature  of  the  development  of  the  internal  plasma  structure  of  a  Hall 
thruster  at  small  timescales.9 

HPHall’s  development  has  been  documented  by  Giuliano  and  Boyd,10  the  results  of  which  are  summarized  here. 
The  HPHall  code  performs  an  axisymmetric  simulation,  commonly  referred  to  as  “hybrid-PIC,”  treating  the 
electrons  via  fluid  approximation  equations  and  treating  the  heavy  species,  namely  xenon  ions  and  atoms,  via  a 
particle-in-cell  (PIC)  method.  The  electron  equations  are  solved  at  a  smaller  time  step,  called  the  electron  subcycle, 
in  order  to  accomodate  the  fact  that  the  speed  of  electrons  is  much  faster  than  that  of  a  xenon  atom  or  ion.  This 
allows  for  electron  fluid  equations  to  be  fully  converged  in  between  xenon  particle  updates.  Hybrid  fluid-particle 
methods  have  been  shown  to  be  successful  in  Hall  thruster  plume  studies1  and  are  computationally  more  efficient 
than  fully  kinetic  methods.11  HPHall,  originally  created  by  Fife  and  Martinez-Sanchez,  was  later  advanced  by 
Gamero-Castano  and  Katz  to  include  such  upgrades  as  a  more  robust  sheath  model  and  a  sputtering  yield  algorithm 
for  the  use  as  an  erosion  model.12  Further  development  and  corrections  were  performed  by  Hofer  et  al.  on  the  heavy 
particle  modeling,  erosion  sub-model,  and  electron  mobility  physics,  continuing  the  development  of  HPHall  to  the 
present  version.8' 13-14  Through  this  past  work,  the  code  has  a  favorable  history  of  presenting  good  agreement  with 
macroscopic  properties,  such  as  discharge  current  and  thrust,  as  well  as  local  properties,  such  as  plasma  density, 
plasma  potential,  and  electron  temperature.  An  example  of  this  time-averaged  output  of  the  code  can  be  seen  in 
Figure  2  which  is  a  representation  of  electron  number  density  at  nominal  operating  conditions,  detailed  in  Table  1. 
For  this  type  of  two-dimensional  contour  plot,  the  anode  is  located  on  the  left  side  and  the  near-field  plume  is 
located  on  the  right  side  of  the  domain. 

HPHall  has  been  developed  to  include  the  capability  to  output  detailed  information  of  each  xenon  species  at 
various  points  in  the  simulation  domain.  This  capability  is  utilized  in  the  current  study  to  improve  the  accuracy  of 
plume  simulation  input  conditions,  specifically  at  the  thruster  exit.  Previous  methods  used  in  determining  these 
conditions  were  semi-empirical.  Comparisons  between  input  conditions  predicted  by  these  two  approaches  will  be 
made  in  Section  IV. 
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Figure  2:  Time-averaged  electron  number  density  contours 


III.2  DSMC-PIC  Methodology 
A.  Overview  and  Collision  Dynamics 

The  simulations  use  a  hybrid  DSMC-PIC  method.  The  DSMC  module  handles  collisions.  The  PIC  module  is 
used  to  move  the  heavier  particles  that  are  influenced  by  the  electric  fields  present.  Finally,  the  electrons  are 
simulated  using  a  detailed  fluid  model. 

The  DSMC  module  handles  collisions  between  heavier  particles  (Xe,  Xe+,  and  Xe++),  such  as  neutral-neutral  and 
ion-neutral  collisions.  The  DSMC  method  uses  virtual  particles  to  simulate  collisions  in  rarefied  gas  flows.  The 
particles  represent  real  ions  and  neutrals  and  are  grouped  in  cells  whose  characteristic  lengths  are  shorter  than  a 
mean  free  path.  Pairs  of  these  particles  are  selected  at  random  and  a  collision  probability  is  evaluated  that  is 
proportional  to  the  product  of  the  relative  velocity  and  collision  cross-section.  This  probability  is  compared  to  a 
random  number  to  determine  if  the  collision  occurs.  If  so,  collision  dynamics  are  performed  to  alter  the  properties 
of  the  colliding  particles.  Two  types  of  collision  dynamics  are  relevant  to  Hall  thruster  plumes:  elastic  (momentum 
exchange)  collisions  and  charge  exchange  collisions  (CEX).  Elastic  collisions  involve  only  exchange  of  momentum 
between  participating  particles.  The  Hall  thruster  plume  is  confined  to  two  different  types  of  momentum  exchange 
collisions,  neutral-neutral  collisions  and  neutral-ion  collisions.  For  neutral-neutral  collisions,  the  variable  hard 
sphere  model  is  employed’.  The  cross-section  for  xenon  is: 

,  ,,  ...  2.12  x  10-18  o  ,  _  „  „ 

oEL{Xe,Xe)  =  — ^ - m  and  a>  =  0.12  (1) 

where  g  is  the  relative  velocity  and  co  is  related  to  the  viscosity  temperature  exponent  for  xenon.  For  neutral-ion 
elastic  interactions,  the  cross-sections  measured  by  Miller  et  al. 1 5  are  used: 


(2) 
(3) 
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t jEL(Xe,Xe+ )  =  (175.26-  27.2  log10(^))  x  1(T20  m 2 
aEL (  Xe, Xe++  )  =  (  103.26  -  17.8  log10 (  g  ))  x  10~20  m2 


Isotropic  scattering  is  assumed  for  both  types  of  elastic  collisions.  Charge  exchange  collisions  pertain  to  the  transfer 
of  one  or  more  electrons  between  an  atom  and  an  ion.  These  cross  sections  are  assumed  to  follow  the  same 
expressions  for  neutral-ion  elastic  collisions.  However,  it  is  also  assumed  there  is  no  transfer  of  momentum 
accompanying  the  charge  exchange,  since  it  is  primarily  a  long-range  interaction. 


B.  Plasma  Dynamics 

The  PIC  module  is  used  to  move  the  heavier  ion  particles  that  are  influenced  by  the  electric  fields,  whereas  the 
lighter  electrons  are  modeled  as  a  fluid.  The  PIC  module  determines  the  charge  density  at  the  nodes  in  the  mesh 
based  on  the  proximity  of  each  particle  to  the  surrounding  nodes.  The  charge  density  is  then  used  to  compute  the 
electric  field  at  each  node.  This  is  accomplished  either  by  incorporating  the  Boltzmann  relation  or  solving  for  the 
potential  directly  using  the  detailed-fluid  model.  The  potential  is  then  differentiated  spatially  to  obtain  the  electric 
fields. 

The  Boltzmann  relation  uses  several  assumptions  applied  to  the  electron  momentum  equation,  such  as  the  fluid 
electron  flow  being  collisionless,  isothermal,  with  no  magnetic  fields  present,  and  that  the  electron  pressure  obeys 
the  ideal  gas  law,  to  arrive  at  the  following10: 


0  (pref  +  efln(nreef)  (4) 

where  </>  is  the  plasma  potential,  <pref  is  a  reference  potential,  k  is  Boltzmann’s  constant,  e  is  the  electron  charge,  Tref  is 
the  constant  electron  reference  temperature,  and  ne  and  nref  are  the  local  electron  number  density  and  reference 
electron  number  density,  respectively.  Reference  values  are  primarily  taken  at  the  thruster  exit  plane.  The 
assumption  of  quasi-neutrality  is  employed  to  obtain  the  electron  number  density  from  the  ion  number  densities. 

Most  of  the  strong  assumptions  made  in  deriving  the  Boltzmann  relation  are  questionable  in  the  plume  under 
consideration,  especially  in  the  near-field.  This  is  due  to  strong  gradients  that  could  make  these  approximations 
inaccurate.  Previous  work  on  Hall  thruster  plume  modeling  shows  the  increased  fidelity  more  detailed  physics 
models  have  over  the  Boltzmann  relation,  which  in  turn  provides  better  agreement  with  experimental  results.2  The 
detailed  electron  model5  increases  the  modeling  fidelity  by  removing  some  of  the  simplifications  made  with  the 
Boltzmann  relation  by  modeling  the  electrons  with  fluid-type  conservation  equations.  The  set  of  Eqs.  (5)-(7)  is  used 
to  increase  the  level  of  accuracy  by  removing  some  restrictions  that  the  Boltzmann  relation  assumes,  e.g.  currentless 
and  isothermal  electrons.  The  relations  are  manipulated  into  useful  forms  for  numerical  simulation  by  introducing  a 
stream-function  (called  the  electron  velocity  potential)  for  the  electron  continuity  equation  and  assuming  steady  state 
for  all  the  equations.  This  results  in  a  set  of  Laplace-type  equations  with  weak  source  terms,  summarized  as  follows: 

Vzi p  —  nenaCi  where  Slip  —  ne  \Te  (5) 
V(<tV0)  =  k-{pSJ2Te  +  oTeV2(ln(ne))  +  aV(ln(ne ))  •  S7Te  +  TeS7o  ■  V(ln(ne))  +  V<r  •  S7Te  )  (6) 

S72Te  —  -V(ln(/ce))  •  S7Te  +  ^-(—j-E+  ^ne(t£  ■  S7)kTe  +  peS7-v^  +  3  ^-veinek(Te  —  Th)  +  nenaCie^j  ^ 

This  system  is  treated  in  the  current  study  in  the  same  manner  as  in  Refs.  17  and  18,  the  results  of  which  are 
summarized  as  follows.  By  treating  the  right-hand  side  as  known  in  Eqs.  (5)-(7),  three  fundamental  plasma 
parameters  are  solved  for,  namely  ve,  </>,  and  Te.  This  is  achieved  by  expressing  the  system  as  a  generalized  Poisson 
equation  and  solving  the  system  with  a  finite  element  solver.  Derivative  calculation  is  handled  by  a  least  squares 
method,  also  presented  in  Ref.  17. 


C.  Boundary  Conditions 

For  computations  of  a  Hall  thruster  plume  in  the  LVTF,  boundary  conditions  must  be  specified  at  the  thruster 
exit,  outflow  surfaces,  chamber  centerline,  and  along  the  thruster  wall,  for  both  the  DSMC-PIC  method  and  for  the 
fluid  electron  model.  Figure  3  shows  a  schematic  of  the  grid  utilized  in  the  current  study.  First,  the  boundary 
conditions  for  the  DSMC-PIC  method  are  presented,  followed  by  those  for  the  fluid  electron  model. 

Some  of  the  macroscopic  properties  of  the  plasma  are  required  at  the  thruster  exit,  namely  the  number  density, 
velocity  components,  and  temperature  of  each  heavy  species  in  the  calculation.  Since  the  particles  exit  the  thruster 
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Figure  3:  Schematic  of  axisymmetric,  unstructured 
grid,  normalized  by  mean  thruster  diameter 


with  an  unknown  radial  velocity  component,  this  results  in  a  velocity  vector  that  is  not  parallel  to  the  center-line  of 
the  chamber.  The  angle  between  the  center-line  and  velocity  vector  is  referred  to  as  the  divergence  angle,  and  since 
it  is  not  known  a  priori,  a  sensitivity  study  is  performed  that  includes  a  range  of  assumed  values.  Particles  that  exit 
the  thruster  plane  are  assigned  a  radial  velocity  component  that  varies  linearly  with  the  distance  from  the  center  of 
the  thruster  channel.  Therefore,  a  particle  exiting  the  center  of  the  channel  will  have  no  radial  velocity  component, 
while  a  particle  exiting  the  channel  near  its  inner  or  outer  wall  will  have  the  highest  radial  velocity  component  and 
the  angle  of  the  velocity  vector  with  respect  to  the  center-line  for  that  particle  will  be  the  divergence  angle  specified. 

The  other  thruster  exit  macroscopic  properties  are  determined  in  general  using  two  methods:  first,  via  a 
combination  of  analysis  and  estimation  in  order  to  match  experimental  operating  conditions  in  the  same  manner  as 
Ref.  2.  The  resulting  DSMC-PIC  input  parameters  produce  thruster  exit  conditions  matching  the  nominal 
experimental  operating  conditions  to  within  5%  (see  Table  1).  This  has  been  noted  in  the  same  work  as  a  source  of 
uncertainty  in  comparison  with  experimental  data.  Therefore,  a  second  method  is  employed:  conditions  are 
obtained  through  analysis  of  HPHall  predictions  at  the  thruster  exit  plane.  This  approach  also  resulted  in  parameters 
matching  experimental  operating  conditions  to  within  5%.  The  particle  properties  predicted  by  each  method  are 
shown  in  Table  2;  see  Section  IV  for  further  discussion  of  these  predictions. 

To  determine  particle  properties  at  the  cathode,  the  assumption  is  made  that  the  mass  flow  consists  solely  of 
neutral  xenon  atoms.  This  assumption  allowed  for  a  clear  calculation  of  boundary  conditions  for  the  fluid  electron 
model,  as  injecting  xenon  ions  from  the  cathode  creates  a  current  density  which  feeds  back  into  determination  of  the 
fluid  boundary  conditions.  The  neutrals  are  assumed  to  have  a  characteristic  temperature  of  1300  K  as  reported  in 
Boyd19.  The  neutrals  injected  at  the  cathode  are  assumed  to  have  sonic  velocity,  which,  using  the  reported  mass 
flow  rate  of  7%  of  the  anode  mass  flow  rate6,  thereby  determines  the  number  density.  The  particle  properties  for 
species  at  the  cathode  are  also  listed  in  Table  2. 

There  is  only  one  solid  surface  that  is  modeled  in  the  present  analysis,  namely  the  thruster  itself.  It  is  assumed  to 
have  a  plasma  potential  of  zero,  i.e.  to  be  electrically  grounded.  All  ions  that  collide  with  the  thruster  wall  are 
neutralized.  For  the  particles  scattered  back  into  the  flow  field  from  the  thruster  wall,  diffuse  reflection  is  assumed 
which  is  characterized  by  the  surface  temperature  of  300  K. 

The  facility  back-pressure  is  modeled  through  static  background  particles.  Each  cell  contains  a  few  particles  with 
velocities  sampled  from  a  zero-centered  Maxwellian  velocity  distribution  function  at  an  assumed  temperature  of  295 
K.  These  particles  participate  in  collisions  with  plume  particles  and  change  the  velocities  of  other  particles,  but  their 
positions  and  velocities  do  not  change.  The  back-pressure  value  for  these  simulations  is  set  to  the  value  near  the 
general  operation  back  pressure, 

1.3  x  10"5  torr. 


6 

46,h  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference 


Thruster  Particle  Properties 

Cathode  Particle  Properties 

Species 

n  (1  /  m3) 

u  (m/  s) 

T  (  K  ) 

n  (1  /  m3) 

u  (m  /  s) 

T  (  K  ) 

Semi- 

empirical 

HPHall 

Semi- 

empirical 

HPHall 

Semi- 

empirical 

HPHall 

Conditions  same  for  both  semi- 
empirical  and  HPHall 

Xe 

1.77  x 
1018 

2.82  x 
1016 

281 

410 

750 

620 

1.78  x 
1020 

320 

1300 

Xe+ 

2.00  x 
1017 

3.78  x 
1017 

19600 

13450 

23200 

45000 

0 

0 

0 

Xe^ 

5.26  x 
1016 

5.62  x 
1016 

27800 

17800 

23200 

105000 

0 

0 

0 

Table  2:  Particle  properties  for  DSMC-PIC 

The  boundary  conditions  for  the  detailed-fluid  model  must  be  determined  in  order  to  solve  Eqs.  (5)  -  (7)  above. 
Since  each  equation  is  Laplace-like,  each  one  requires  specification  of  either  a  Dirichlet  (direct)  value,  or  a  von 
Neumann  (gradient)  value.  The  potential,  (f),  the  electron  temperature,  T„  and  the  quantity  nrxf,  are  specified  as 
either  direct  or  gradient  values  at  each  boundary  in  the  simulation.  The  value  to  which  each  is  set  is  specified  in 
Table  3.  Each  boundary  condition  is  direct  except  for  the  following:  1)  each  boundary  condition  for  the  axis  of 
symmetry  is  a  gradient-type  condition,  2)  the  wall  of  the  thruster  has  a  gradient-type  condition  for  the  electron 
temperature,  and  3)  the  outflow  regions  have  a  gradient-type  condition  for  neve.  Due  to  the  grounding  of  the 
thruster,  the  plasma  potential  condition  is  direct  and  set  to  zero  for  the  thruster  walls  and  set  to  reference  values  for 
the  thruster  exit  plane  and  cathode.  The  electron  temperature  is  set  to  characteristic  values  based  on  measurements 
in  Ref.  6.  The  boundary  conditions  for  neve  are  determined  using  conservation  of  current  as  follows: 

Id  ~  4  +  h 

Ij  —  ]\  *  Ate 

h 

Jd,  cathode  a 

cathode 

where  (8)  and  (9)  are  used  to  compute  current  density  at  the  thruster  exit  and  (10)  computes  current  density  at  the 
cathode,  from  which  neve  is  found  by  dividing  by  the  elementary  charge  e.  The  sign  convention  results  from 
electrons  flowing  out  of  the  cathode,  towards  the  anode.  Therefore,  the  electron  velocity  stream  function  must  be 
positive  at  the  cathode  and  negative  at  the  thruster  exit. 


(8) 

(9) 

(10) 


Boundary 

Condition 

Variable 

Thruster 

Cathode 

Symmetry-line 
(all  gradient-type) 

Outflow 

Exit  Plane 

Body 

dW 

-4.64  x  1022  m"2  s'1 

0  m'2  s"1 

2.85  x  1023  nf2  s'1 

0  mV 

0  m"3  s"1 
(gradient) 

<t> 

40  V 

0  V 

12  V 

OVm'1 

1  V 

T 

±  e 

25  eV 

0  eV  nf1 
(gradient) 

15  eV 

0  eV  m"1 

1  eV 

Table  3:  Detailed  fluid-electron  boundary  conditions 


7 

46th  A1AAJASME/SAE/ASEE  Joint  Propulsion  Conference 


IV.  Results 

The  HPHall  simulations  are  performed  as  described  in  Ref.  10:  neutral  xenon  atoms  are  injected  for  20,000  time 
steps,  at  which  point  xenon  ions  are  also  injected  for  an  additional  80,000  time  steps.  The  time  step  size  used  for  the 
xenon  particles  is  5  x  10  s  seconds,  with  an  electron  subcycle  time  step  of  5  x  10’11  seconds.  The  HPHall 
simulations  used  approximately  300,000  particles.  The  plume  simulations  presented  use  a  total  of  500,000  particles 
at  steady  state  over  a  domain  of  3,925  triangular  cells.  The  simulation  runs  for  30,000  time-steps  to  reach  a  steady 
state  and  then  for  another  10,000  time-steps  to  sample  macroscopic  data.  The  time-step  size  is  1  x  10"5  seconds, 
resulting  in  a  total  sampling  time  of  1  second.  Cases  are  run  at  a  maximum  divergence  angle  of  10°. 

A.l  HPHall  and  Semi-empirical  Predictions  of  Thruster  Exit 

Shown  in  Figures  4-6  are  velocity  distribution  functions  (VDF)  based  on  the  semi-empirical  method  and  from 
HPHall  for  each  xenon  species.  These  conditions  are  then  used  to  construct  a  Maxwellian  VDF  which  in  turn 
generates  injected  particles.  The  semi-empirical  method  predicts  these  boundary  conditions  through  estimation 
based  on  thruster  operating  conditions.  The  method  also  requires  certain  assumptions  be  made,  such  as  specifying  a 
priori  the  fraction  of  xenon  atoms  ionized.  The  main  disadvantage  of  the  semi -empirical  method  is  its  imprecision: 
assumptions  must  be  made  about  flow  conditions  with  no  knowledge  regarding  the  physical  situation.  This 
disadvantage  was  noted  in  previous  work.2  As  an  alternative  approach,  using  HPHall,  the  xenon  species  are 
sampled  along  the  thruster  exit,  with  results  shown  alongside  the  semi-empirical  results  in  Figures  4-6.  The  VDF 
for  each  species  predicted  by  HPHall  is  structurally  different  in  all  cases;  however,  each  VDF  predicted  by  HPHall 
closely  matches  the  structure  of  a  Maxwellian  VDF.  This  congruence  allows  for  the  calculation  of  new  boundary 
conditions  for  the  plume  simulations  while  maintaining  the  assumption  of  a  Maxwellian  VDF  at  the  thruster  exit. 
The  new  boundary  conditions  are  determined  by  constructing  Maxwellian  VDF’s  that  align  with  those  predicted  by 
HPHall;  see  Figures  7-9. 


Cj,  m/s 

Figure  4:  Semi-empirical  prediction  (solid  line) 
and  HPHall  prediction  (dashed  line)  for  xenon 
neutrals 


Cj,  m/s 

Figure  5:  Semi-empirical  prediction  (solid  line) 
and  HPHall  prediction  (dashed  line)  for  single - 
charged  xenon  ions 
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OP  (tD)J  t3P  (t3)j 


0  10000  20000  30000  40000 

C1#  m/s 

Figure  6:  Semi-empirical  prediction  (solid  line) 
and  HPHall  prediction  (dashed  line)  for  double- 
charged  xenon  ions 


Cj,  m/s 

Figure  7:  Constructed  Maxwellian  (solid  line) 
and  FIPFIall  prediction  (dashed  line)  for  xenon 
neutrals 


0  10000  20000  30000  40000 

Cj,  m/s 

Figure  8:  Constructed  Maxwellian  (solid  line) 
and  HPHall  prediction  (dashed  line)  for  single- 
charged  xenon  ions 


0  10000  20000  30000  40000 

Clf  m/s 

Figure  9:  Constructed  Maxwellian  (solid  line)  and 
HPHall  prediction  (dashed  line)  for  double- 
charged  xenon  ions 
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A  comparison  of  mean  velocity,  temperature,  and  number  density  between  the  semi-empirical  method  and  the 
constructed  Maxwellian  is  made  in  Table  2  above.  When  compared  to  the  VDF’s  constructed  through  HPHall,  the 
semi-empirical  method  generally  over-predicts  the  mean  velocity  of  ions  exiting  the  thruster,  while  under-predicting 
the  mean  velocity  of  neutrals.  The  xenon  ion  VDF’s  predicted  by  the  semi-empirical  method  are  also  narrower  than 
the  constructed  VDF’s:  this  results  in  the  FlPHall-constructed  input  conditions  specifying  a  much  higher 
temperature  than  the  temperature  assumed  by  the  semi-empirical  method.  Note  that  the  trend  is  reversed  for  the 
neutral  xenon  species:  the  constructed  VDF  is  narrower  than  the  semi-empirical  method,  resulting  in  a  lower 
specified  temperature  for  the  neutrals. 

A.2  Plume  Simulations 

There  is  a  wide  disparity  between  the  two  input  conditions,  and  as  a  result,  the  plume  simulations  that  utilized 
each  set  of  conditions  were  also  disparate.  Figures  10  and  1 1  illustrate  the  effect  of  the  different  thruster  conditions 
on  general  features  of  the  plume  simulations.  The  differences  in  the  current  density  contours  are  mainly  attributed  to 
the  predicted  axial  velocity  for  each  method.  While  the  mass  flux  in  each  case  is  the  same,  the  semi-empirical 
method  predicts  an  axial  velocity  greater  than  that  predicted  by  the  constructed  method  by  a  factor  of  about  1.5  for 
both  single  and  double  charged  xenon  species.  Because  the  reference  potential  is  kept  constant  in  these 
comparisons,  the  difference  in  axial  velocity  is  primarily  responsible  for  the  longer  beam-like  plume  for  the  semi- 
empirical  conditions,  as  seen  in  Figure  10.  The  effect  of  these  changes  on  the  electron  temperature  can  be  seen  in 
Figure  1 1 :  the  constructed  conditions  result  in  a  higher  peak  electron  temperature  near  the  thruster  exit. 


Z  (  D  1 

Figure  10:  Current  density  comparison, 

constructed  input  condition  (top)  and  semi- 
empirical  input  condition  (bottom) 


Z  (  D  ) 

Figure  11:  Electron  temperature  comparison, 
constructed  input  condition  (top)  and  semi- 
empirical  input  condition  (bottom) 


B.  Experimental  Comparison 

Measured  data  are  reported  in  Ref.  6,  where  data  acquisition  methods  are  also  described.  All  measurements  are 
taken  at  the  nominal  thruster  operating  condition.  Model  assessment  is  performed  through  comparison  with  data 
measured  by  a  nude  Faraday  cup  probe.  For  the  operating  condition  under  consideration,  over  64,000  individual  ion 
current  density  measurements  were  taken  in  the  near-field  region.  The  probe  did  not  incorporate  a  guard  ring  in  this 
case,  as  the  local  Debye  length  made  this  impractical.  Measurement  uncertainty  is  well  documented  in  Ref.  6. 
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Therein  uncertainties  due  to  probe  collection  area,  sheath  expansion,  and  cosine  losses  are  reported  and  conservative 
estimates  suggest  uncertainties  of  ±10%  for  ion  current  density  measurements. 

Shown  below  in  Figures  12  and  13  are  contour  plots  comparing  experimental  measurements  and  various 
simulations.  Contour  plots  are  used  to  show  the  effect  of  method  used  to  determine  thruster  exit  conditions.  The 
main  structural  similarities  between  experimental  measurements  and  simulations  are  in  and  around  the  highly 
collimated  plume  at  the  thruster  exit.  Throughout  the  first  thruster  diameter,  the  main  thruster  plume  remains  beam¬ 
like  regardless  of  which  method  is  used  to  determine  thruster  boundary  conditions.  Each  method  also  predicts 
coalescence  of  the  main  beam  near  the  chamber  centerline.  However,  the  semi-empirical  method  over-predicts  the 
peak  current  density  in  the  main  ion  beam  by  33.79%,  whereas  the  constructed  boundary  conditions  come  much 
closer  to  this  peak,  differing  by  only  -1.26%.  This  discrepancy  is  highlighted  in  Figure  14,  which  shows  current 
density  variation  along  the  discharge  chamber  centerline.  Figure  14  illustrates  the  advantage  of  using  HPHall  to 
construct  boundary  conditions  that  are  a  closer  approximation  to  the  physical  state  at  the  thruster  exit  than  the 
approximation  using  the  semi-empirical  method. 

Figures  15-18  illustrate  differences  in  plume  structure  by  examining  radial  variation  of  ion  current  density  at 
different  locations  downstream  of  the  thruster  exit.  Figures  15  and  16  compare  current  density  measurements  and 
predictions  via  the  semi-empirical  method,  whereas  Figures  17  and  18  compare  the  same  measurements  with 
predictions  via  the  constructed  method.  There  are  two  peaks  in  current  density  immediately  downstream  of  the 
thruster  exit  (0.1  D),  due  to  the  mass  flows  from  both  the  thruster  and  the  cathode.  Both  simulations  and 
measurements  show  that  further  downstream  the  two  peaks  begin  to  coalesce,  with  measurements  having  almost  no 
sign  of  a  second  peak  by  1.2  D  downstream  of  the  thruster  exit.  Both  simulations  predict  this  coalescence  occurring 
in  a  shorter  downstream  distance  than  measurements:  by  around  0.5  D  downstream  of  the  thruster  exit,  both 
simulations  predict  a  single  peak  current  density. 


0  0.5  1  1.5  2  2.5  3 


Z  (  D  ) 


0  0.5  1  1.5  2  2.5  3 


Z  (  D  ) 


Figure  12:  Current  density  comparison, 

measured  data  (top)  and  semi-empirical  input 
condition  (bottom) 


Figure  13:  Current  density  comparison, 

measured  data  (top)  and  constructed  input 
condition  (bottom) 
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J  (  mA / cm  2 


Figure  14:  Current  density  variations  along  the  discharge 
chamber  centerline,  measured  from  the  thruster  exit  (note 
error  bars  on  measured  data  are  ±10%) 


Figure  15:  Radial  current  density  comparisons 
between  measurements  and  predictions  via  semi- 
empirical  method  at  various  axial  stations 


Figure  16:  Radial  current  density  comparisons 
between  measurements  and  predictions  via  semi- 
empirical  method  at  different  axial  stations 
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R  (  D  ) 


Figure  17:  Radial  current  density  comparisons 
between  measurements  and  predictions  via 
constructed  method  at  various  axial  stations 


R  (  D  ) 


Figure  18:  Radial  current  density  comparisons 
between  measurements  and  predictions  via 
constructed  method  at  different  axial  stations 


Neither  simulation  accurately  predicts  the  plume  structure  outside  of  the  near-field  plume  region,  as  seen  in 
Figures  15-18.  This  is  due  to  two  major  contributing  factors.  First,  the  numerical  model  is  limited  in  requiring 
fluid-type  boundary  conditions  throughout  the  entire  domain,  which  are  typically  not  available.  For  example,  Eq. 
(7)  requires  a  direct  value  for  electron  temperature  set  at  the  outflow  surfaces,  but  acquisition  of  electron 
temperature  measurements  for  such  a  large  range  of  spatial  locations  is  not  feasible.  Second,  the  current  numerical 
model  makes  strong  assumptions  about  the  makeup  of  the  cathode  mass  flow.  Currently,  only  xenon  neutrals  are 
injected  from  the  cathode.  On  comparing  the  computational  predicted  current  density  at  the  cathode  with  the 
measured  current  density,  it  is  clear  that  the  cathode  flow  also  contains  some  xenon  ions,  as  the  measured  current 
density  is  nearly  twice  that  of  the  predicted  current  density.  Previous  work  on  cathode  modeling  also  indicates  the 
presence  of  xenon  ions19.  The  current  study  must  assume  a  cathode  mass  flow  consisting  of  only  neutrals  because 
there  is  no  data  available  on  the  ion  properties  at  the  cathode  outside  of  current  density  measurements. 


V.  Conclusions 

A  general  purpose,  hybrid  DSMC-PIC  code  has  been  supplemented  with  a  simulation  package,  HPFlall. 
HPFlairs  capability  to  predict  thruster  exit  conditions  has  been  assessed.  The  refined  conditions  more  accurately 
characterize  what  was  previously  an  area  of  high  uncertainty,  as  determined  in  Ref.  2.  The  new  approach  provides  a 
scientific  alternative  to  the  semi -empirical  estimation  of  particle  properties  at  the  thruster  exit  by  analysis  of  velocity 
distribution  functions.  These  refined  input  conditions  for  plume  simulations  enhance  simulation  accuracy  in  the 
near-field  plume  region,  predicting  peak  current  density  at  the  thruster  exit  to  within  1.25%. 

One  area  of  future  study  to  better  understand  the  physical  situation  in  the  plume  is  the  characterization  of  the 
xenon  species  at  the  cathode.  As  has  been  previously  noted,  the  mass  flow  at  the  cathode  is  made  up  of  xenon 
neutrals  and  ions.  Since  there  is  a  large  flux  of  ions  at  the  cathode  that  is  not  currently  being  modeled,  further 
investigation  into  the  nature  of  the  cathode  flow  is  required  to  improve  simulation  accuracy  outside  the  near-field 
region  of  the  plume. 
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